Reconstructing nonlinearities with intermodulation spectroscopy 
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We describe a method of analysis which allows for reconstructing the nonlinear disturbance of a 
high Q harmonic oscillator. When the oscillator is driven with two or more frequencies, the nonlin- 
earity causes intermodulation of the drives, resulting in a complicated spectral response. Analysis of 
this spectrum allows one to approximate the nonlinearity. The method, which is generally applicable 
to measurements based on resonant detection, increases the information content of the measurement 
without requiring large detection bandwidth, and optimally uses the enhanced sensitivity near res- 
onance to extract information and minimize error due to detector noise. 



A harmonic oscillator with high quality factor Q, i.e., 
large stored energy in relation to the energy lost in each 
cycle, is a very useful tool for precision measurement. 
When the oscillator interacts with an object of interest, 
the first order correction to the linear response is a shift 
of the resonant frequency cjo, which results in a large 
change of the amplitude and phase of the response near 
resonance. Resonant detection is employed in a wide va- 
riety of measurements, and particularly in nanoscience, 
where examples include the radio frequency single elec- 
tron transistor [l[ , dispersive readout of superconducting 
qubits [2], measuring the deflection of a nanomechani- 
cal oscillator [3], and nanoscale imaging of surfaces with 
dynamic atomic force microscopy (AFM) [4| . In this let- 
ter we describe a method of resonant detection which 
exploits the nonlinear phenomenon of intermodulation, 
to generate a spectral response near resonance that en- 
hances the information content of the measurement. We 
provide a general theoretical framework for analyzing 
the intermodulation spectrum and for reconstructing the 
nonlinearity. 

The sensitivity of resonant detection can be under- 
stood from the large transfer gain, G(cj), of a high-Q 
oscillator. In the narrow frequency band B = uoq/Q cen- 
tered at resonance, the oscillator gives large response to a 
small stimulus. When the linear oscillator interacts with 
the object of interest, the equation of motion becomes 
nonlinear and higher harmonics of the drive frequency 
will appear in the response. One could in principle an- 
alyze the amplitude and phase of these higher harmon- 
ics to extract the nonlinearity. The problem with this 
approach is: (i) it requires a very large detection band- 
width, £?d = nido to get n harmonics and (ii) the trans- 
fer gain for harmonics is far less than one, falling off as 
1/uj 2 above resonance, so the response at higher harmon- 
ics is often buried in the detector noise. In the field of 
AFM, the latter problem was addressed by coupling more 
oscillators, utilizing additional torsional {jjj and flexural 
[6|,l3] eigenmodes of the cantilever. In this context, differ- 
ent types of multi-tone excitation schemes have been put 
forward as a means of extracting qualitative information 
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FIG. 1: (color online). Intermodulation spectroscopy: two 
pure tones drive the resonator which interacts with an ob- 
ject of interest. The detected response spectrum consisting 
of many intermodulation products can be inverted to recon- 
struct the nonlinear interaction. 

about the nonlinearity 

An alternative approach, sketched in Fig. [TJ is to drive 
the oscillator with two pure tones uj\ and so as to pro- 
duce intermodulation products (IMPs) or "mixing prod- 
ucts" which occur at the frequencies nui + mu2 (n and 
m integers), where the order of the IMP is defined by 
| n | + \m\. The drive frequencies can be chosen in a va- 
riety of ways, with the objective of creating many IMPs 
near resonance where the sensitivity is enhanced. These 
IMPs measured near resonance represent a partial spec- 
tral response of the system from which one would like to 
reconstruct the nonlinearity. 

Our work is motivated by the desire to understand 
the information content of intermodulation spectra, as 
measured in superconducting microresonators at GHz 
frequencies [ll|, [l2j , and in AFM cantilevers oscillating 
at several hundred kHz while interacting with a surface 
[13|, LL4( (see Fig. [2]). However, we wish to stress that the 
reconstruction algorithm described here can be applied 
to a wide variety of experiments which exploit resonant 
detection. To test the reconstruction algorithm, we use 
simulated data obtained by numerical integration of the 
equation of motion of the nonlinear oscillator, where, in 
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FIG. 2: (color online). Measured intermodulation spectra 
taken from a nonlinear superconducting resonator [l2|] (left 
panel) and an AFM cantilever interacting with a surface [3] 
(right panel) . In these experiments the two central peaks with 
largest response are the drive tones. The amplitude in dB is 
normalized to the maximum response. 

contrast with experiment, we can work backward from 
the simulated intermodulation spectrum to the known 
expansion coefficients of the nonlinearity, thus allowing a 
more systematic study of errors. 

The equation of motion of a single-mode harmonic os- 
cillator, which is driven by two pure tones, subject to 
viscous damping, and perturbed by a nonlinearity, reads: 

C + ^ + C = ^m(C) + T x cos(n lT ) + T 2 cos(n 2 r) . (1) 

Here, Qi j2 = ^1,2/^0? dots denote derivatives with re- 
spect to the dimensionless time r = uo^t^ and the di- 
mensionless coordinate ( is the deviation of the oscilla- 
tor from the equilibrium position, suitably normalized. 
The dimensionless nonlinear "force" .Fni(C) is assumed 
conservative, with no explicit time-dependence. 

While the motion of a damped, driven nonlinear sys- 
tem of the type $1]) can in general be either regular or 
chaotic, we restrict ourselves here to nonlinear it ies and 
drive strengths which are weak enough, such that for 
commensurate drive frequencies with greatest common 
divisor = gcd(^i, ^2) the motion in the steady state 
is periodic in 2ir/AQ. We may therefore expand the re- 
sponse in a discrete Fourier series, £(r) = J2 k Ck^ zk AQ 7 \ 
where the expansion coefficients Ck are the complex num- 
bers (amplitude and phase) which constitute the inter- 
modulation spectrum. With Fourier expansions of the 
nonlinear force J r n i(C( r )) = ^2k^ nl ^ e%k AQ 1 "> an ^ ^he 
drive force, J^di^) = ^2k^~ d ^ elk A ^ T > ^ ne equation of 
motion (pQ) becomes 

Ck = Gk{Fd,k + ^nl,/c) = Cfc^ + Gk^rA.k , (2) 

where G k = 1/(1 - k 2 AQ 2 + ikAQ/Q). Note that while 
([2]) is an exact equation, it is not a solution of (pQ), since 
^ni,/c depends on all the Ck in a highly nonlinear way. 

To proceed, we represent the nonlinear force as a poly- 
nomial in C(t), Jni(CW) = lL,JLi9j( J ( r ), which allows 



us to express the Fourier components of the nonlinear 
force ^ni^ as a linear combination of the polynomial co- 
efficients, gy. 

00 

^ni,fc = ^2/H kj gj, (3) 

where the matrix element H k j is the kih spectral com- 
ponent of the Fourier transform of C J ( r )? alternatively 
written as H kl = Ck and H k j+i = J2 k > Ck'Hk-k'j- 

For a known nonlinear force, we immediately find a 
first approximation of the spectrum, by evaluating 
^ni,/c in © at the free solution = J^/c- In the next 
step we use in the same way to obtain Ci?\ etc. Each 
iteration generates new spectral response, but only at the 
intermodulation frequencies. If this procedure converges 
to a solution Ck: then our assumption of periodic orbits is 
certainly justified. However, here we are interested in the 
inverse problem of finding the nonlinear force, starting 
from Cfc XP ? which is measured with limited bandwidth in 
the presence of noise. Inspection of Eqs. (|3]) and ([2j) 
reveals that one must invert the matrix H to find the 
expansion coefficients gj of the nonlinearity: 

^•^(H- 1 )^^-^ for jG{l,..., Jmax }. (4) 
k Gk 

The inverse problem is generally complicated by the 
fact that the measured spectrum £jl xp has only a finite 
number N p of peaks that can be observed above the de- 
tector noise level, notably those close to the resonance 
frequency. Nevertheless, with a proper intermodulation 
drive scheme one can achieve enough intermodulation 
peaks with good signal-to-noise ratio to construct a ma- 
trix H, and calculate the coefficients gj. We use the 
following algorithm: 

1. Choose the drive frequencies so that the frequency 
spacing of IMPs, A^, gives N p peaks in the measured 
intermodulation spectrum ^ xp at frequencies kAQ for 
k G {fci, . . . , fcjVp}? where N p is at least twice the number 
of expansion coefficients j max desired in approximating 
the nonlinearity. 

2. Define the intermodulation response Ck to be Ck = (| xp 
for k G {fci, . . . , fc/Vp} and Ck = otherwise. 

3. With the intermodulation response Ck: calculate the 
matrix H, restricted to j < j max and k < /c max • j max , 
where /c max = max{fei, ...,k Np }. 

4. Calculate the j m ax expansion coefficients gj (re- 
stricted to be real) using Eq. (|4|), with the inverse re- 
placed by a pseudo inverse, the sum restricted to k G 
{fci, . . . , fciVp}? and H restricted to k G {k\, . . . , kjsr p } 
and to j < jmax- The pseudo inverse gives a least 
square fit of the j max coefficients gj to the N p equations 

Ck = Cfc 0) + Gk YsflT Hkjgj, thus reducing sensitivity to 
both measurement errors and systematic errors in the 
calculation of H. 
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We have tested this reconstruction algorithm using re- 
alistic parameters for two different experimental realiza- 
tions. In each case, the response of the nonlinear oscilla- 
tor was simulated by numerical integration of Eq. (pQ) us- 
ing the Fortran solver DDASKR [15]. The output of this 
integrator was sampled appropriately and fast Fourier 
transformed to generate an intermodulation spectrum. 
Tests were made to ensure that there was no Fourier leak- 
age of the spectral peaks, so that the background level in 
the spectrum corresponded to the error tolerance set in 
the integrator. Detector noise was simulated by adding 
random gaussian noise A£(r) to each time sample, before 
the Fourier transform. 

As a first example we consider a superconducting 
coplanar wave-guide resonator with a Josephson tunnel 
junction or weak link, currently studied for low- noise am- 
plification and the readout of quantum bits Jj], fl6l - [l9| . 
The classical dynamics of the resonator can be mapped 
to our model equation by using £ = I/Iq, where / is the 
current in the center strip and io the critical current. 
The quality factor Q can be engineered in a wide range, 
up to 10 6 . Typically a nonlinear kinetic inductance pro- 
vides a nonlinear "force" which is constrained by gj = 
for even j, because the inductance dF^x/dQ should not 
depend on the sign of the current (. We can thus employ 
a drive scheme which is only sensitive to IMPs of odd or- 
der, with drive frequencies closely spaced near resonance, 
Aft Efi 2 -Oi <Cl. This drive scheme was used for the 
experimental spectra shown in Fig. 2. For weak drive, 
spectra are well described by the lowest nonlinearity of 
Kerr- type, #3, but for strong drive, higher nonlinear co- 
efficients (#5 ,#7, ...) play a role. 

In Tab. U we give the percent error in the coefficients 
grecon w j 1 i c j 1 were reconstructed from data that was simu- 
lated using three coefficients #3, #5, gj. Results are shown 
for two different simulated noise levels (average of 1000 
simulations) and two different quality factors. In the ab- 
sence of noise, AC = 0, the small errors in the reconstruc- 
tion are explained by the fact that we only used 20 closely 
spaced spectral peaks around the resonance, leading to 
a slight miscalculation of the matrix H. These system- 
atic errors, which are typically small, can be reduced by 
using an improved reconstruction algorithm which esti- 
mates the unmeasured spectral peaks to achieve a self- 
consistent solution of Eq.(|2|). In the presence of noise, 
reconstruction from the intermodulation spectrum gives 
excellent results for all coefficients. For comparison, we 
applied the reconstruction algorithm to a simulated spec- 
trum of the same number of harmonics, generated with a 
single drive having strength J 7 ^ = J r di+J r d2? fc> r the same 
parameters given in Tab. HI In the presence of noise, the 
spectrum of harmonics was unable to reconstruct the co- 
efficients #5 and #7. This nicely illustrates the advantage 
of resonant detection in an intermodulation scheme. 

As a second example, we consider an application which 
contains coefficients gj of both odd and even order. In 
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TABLE I: Relative errors Ag 3 = \( 9j - gf cori )/gj\ in %, for 
different strengths of random Gaussian noise in time, A£(r), 
relative to the maximum response amplitude. Entries 0.00 
mark relative errors smaller than 5 • 10 -5 . To simulate data, 
we used #1 = 0,g 3 = 10~ 3 , g 5 -10~ 4 , g 7 10~ 5 . These 
four coefficients were free parameters for the reconstruction. 
The two drives were centered on resonance. For the case 
Q = 50, AQ = 1/499.5 and F dl = ^2 = 0.03. For the case 
Q = 500, AQ = 1/4999.5 and F dl = ^2 = 0.003. In both 
cases, 2 • 10 5 time samples were used in the Fourier transform. 



dynamic atomic force microscopy (AFM), the resonator 
is a cantilever oscillating about its equilibrium position 
with dimensionless amplitude C = z/z s , effective mass m, 
and stiffness k = muj^. The nonlinearity is provided by 
the tip-surface force, F n \ = F ts (C), which we wish to de- 
termine. Intermodulation spectroscopy can be performed 
while scanning over a surface, measuring the dominant 
IMPs at each image point [HI, Q3] • The inversion algo- 
rithm described here can be used to rapidly determine the 
force-distance curve at each image pixel, while scanning 
at normal speeds for dynamic AFM, one of the major 
objectives of AFM development [il. l20|. 

The nonlinear force typically encountered in AFM ex- 
periments has an attractive region close to a repulsive 
surface. We restrict our study to conservative forces, 
derivable from a potential function, F ts = —dU/dz, and 
desire to reconstruct a polynomial approximation of the 
dimensionless force J-ts(C) — F ts (z + Zo)/(z s h), where 
zq is the probe height above the surface in the absence 
of drive signals, an experimentally tunable parameter. 
For a Morse potential, this force has the functional form 
J^ts = 2a(e- 2 ^- a ^ x - e-^- a ^ x ), which is described by 
the location of the potential minimum, a, its width, A, 
and an overall scale factor, a [21]. We used this nonlin- 
ear force to simulate data with a quality factor Q = 50 
and two drives centered around resonance. Because this 
drive scheme generates near resonance only IMPs of odd 
order, intermodulation peaks around twice the resonance 
frequency were needed to gain information about the co- 
efficients gj with even j. 20 peaks around resonance and 
10 peaks around 2uo were taken for the reconstruction. 
All other spectral information was set to zero. 

The reconstruction algorithm was run using the ampli- 
tude and phase of these N p = 30 peaks, and j max = 15 
coefficients gj were taken for the polynomial represen- 
tation of the nonlinear force. The result of the recon- 
struction is shown in Fig. [3l where we see excellent cor- 
respondence between the original Morse force curve and 
the reconstructed polynomial. 

To make the inversion more realistic, we added 
Gaussian noise corresponding to a maximum-spectral- 
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FIG. 3: (color online). The original force i 7 ts(C) used in the 
simulation, a plot of the polynomial described by the recon- 
structed coefficients gj, and the latter in the presence of noise. 
Inset: The spectrum with noise (blue) and overlaid without 
noise (red). Parameters: a = 0.01, a — —0.6, A = 0.1, 
Aft = 1/499.5, and F dl = T d2 = 0.007. 

amplitude-to- noise ratio ~ 10 -4 (Fig. 3, inset), a dy- 
namic range typical of the optical lever detectors in most 
AFMs. With this noise we see that the polynomial miss- 
represents the force at maximum distance from the sur- 
face, where small changes in the coefficients gj with large 
j have a large effect. We also note that the polynomial 
may in general miss-represent the force near the point of 
maximum curvature, as many coefficients are needed to 
reconstruct sharp kinks in the force-distance curve. To 
improve upon these results, we can carry the reconstruc- 
tion one step further. Using prior knowledge of the inter- 
action in a particular experiment, for example that F ts 
goes asymptotically to zero for large £ and that a sharp 
kink exists at the contact point, one can make a weighted 
fit of the reconstructed polynomial to a particular force 
model containing far fewer than j max parameters. 

Thus it is possible to extract a very nonlinear force- 
distance curve in the fast scanning mode without using 
high frequency components of the cantilever response 
spectrum. Intermodulation spectroscopy allows one to 
collect information in the form of coefficients gj with 
large j, by measuring high order IMPs at much lower fre- 
quencies than the j th harmonic, where such information 
would occur with a single drive tone. The ability of two 
drives to down-convert this information to lower frequen- 
cies, where it can be put near resonance and acquired 
with good transfer gain, is the major strength of the in- 
termodulation spectroscopy technique. The method will 
be accurate if the dominant part of the spectral infor- 
mation is contained in the peaks which are used to con- 
struct the matrix H, and if the traversed force curve is 
well approximated by j max coefficients. Note that both 
amplitude and phase data are necessary to calculate the 
matrix H. 



We briefly note that we also successfully performed 
a simulated reconstruction of force-distance curves with 
two other drive schemes: One with a low- frequency drive 
and a drive at resonance, as sketched in Fig. 1, and 
another scheme with fti = 1/3 and ft2 = 2/3 + Aft. 
These drive schemes have the advantage that they pro- 
duce IMPs of both odd and even order near resonance, 
and are therefore better with the highest Q oscillators, 
and for nonlinearities containing both odd and even co- 
efficients. 

In summary, intermodulation spectroscopy is a mea- 
surement technique which employs a high Q oscillator 
and two drive frequencies to extract response caused by 
a nonlinear disturbance. The drive frequencies are cho- 
sen commensurate with a greatest common divisor Aft, 
and placed so as to produce many IMPs near resonance, 
where large transfer gain allows for detection with good 
signal-to- noise ratio. A reconstruction algorithm was de- 
rived and tested against simulated data including noise. 
Excellent reconstruction of the nonlinearity was demon- 
strated with relatively few spectral peaks obtained at low 
frequencies. Together with the reconstruction algorithm, 
intermodulation spectroscopy should find wide-ranging 
use with the many types of sensors that are based on 
disturbance of a resonance. 
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dation for Strategic Research. 
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